function result = gauss_int1d(function_f, a, b, N_gauss_int1d)

[gauss_points, gauss_weights] = gauss_points_weights_1d(N_gauss_int1d);

A = 0.5*(b-a);
B = 0.5*(a+b);

sum = 0.0;
if mod(N_gauss_int1d,2) == 0
    for i = 1:N_gauss_int1d/2
        sum = sum + gauss_weights(i)*(function_f(B+A*gauss_points(i)) + function_f(B-A*gauss_points(i)));
    end
    result = A*sum;
else
    sum = sum + gauss_weights(1)*function_f(B);
    for i = 1:(N_gauss_int1d-1)/2
        sum = sum + gauss_weights(i+1)*(function_f(B+A*gauss_points(i+1)) + function_f(B-A*gauss_points(i+1)));
    end
    result = A*sum;
end

end